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We have developed a new tool for numerical work in General Relativity: GR- 
workbench. While past tools have been ad hoc, GRworkbench closely follows the 
framework of Differential Geometry to provide a robust and general way of com- 
puting on analytically defined space-times. We discuss the relationship between 
Differential Geometry and CH — h classes in GRworkbench, and demonstrate their 
utility. 

1 Introduction 

We have developed a new class of computational tool for General Relativity. Pre- 
vious tools have fallen into three categories; large scale simulations that evolve 
space-times from initial conditions, symbolic manipulators, and ad hoc numerical 
systems. p. 

GRworkbenchd uses numerical variants of standard differential geometric enti- 
ties to rigorously define space-times in a fashion amenable to computation. This 
system forms a strong base on which to build generally applicable numerical algo- 
rithms, capable of acting on any space-time for which a basic analytic definition is 
available. 

This paper will focus on GRworkbcnch's roots in differential geometry and will 
demonstrate the software's wide applicability to problems via consideration of a 
specific example — geodesic tracing in the Schwarzschild space-time. For a discussion 
of numerical algorithms, visualization teiJmiques and the user interface employed 



2 Discrete differential geometric structure 

We follow the conventions of Hawking and Ellis§: A C n- dimensional manifold 
is a set A4 together with a C atlas {{Ua,4>a)}, that is to say a collection of 
charts {U^, 4>a) where the Ua are subsets of M. and the 4>a are one-to-one maps of 
the corresponding Ua to open sets in R" such that 

1. the Ua cover M, i.e. M = Uq^o; 

2. if Ua n Up is non-empty, then the map 

is a C map of an open subset of R" to an open subset of R". 

It is assumed, as in Hawking and Ellis, that we are dealing with paracompact, 
connected, C°° HausdorfJ manifolds without boundary. 
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Figure 1. A precessing orbit near the event horizon of a Schwarzschild black hole (M = 0.5), 
computed and visualized in GRworkbcnch. The fiat protrusion from the event horizon is the 
(/< = 0, 27r boundary of the spherical polar chart. 



A computer cannot numerically represent the set A4 — that is, it cannot repre- 
sent A4 using numbers. Nor can it so represent Ua or the mapping 0^ for any chart 
{bla,4'a)- The manifold, subsets of the manifold, and functions on the manifold 
are abstract entities; the computer cannot deal with them numerically. This is not 
to say that computers cannot deal with these abstract entities symbolicaUy — there 
exist numerous symbolic manipulators, such as MathematicaB and SheepEl, for this 
purpose. 

Computers can, however, operate numerically in R". The set 4>a{l^a) can be 
represented numerically, as can the function (pa o (pj^^ ■ (f>p{Uar]U/3) — > (fiaiUa <^l^p)- 
Computers represent real numbers as floating-point numbers. This is analogous to 
base-2 scientific notation, where (a, b) represents a x 2^. On a modern computer, 
real numbers are typically represented using 64 bits, meaning that the computer can 
represent up to 2^^ sa 10^^ rational numbers, spaced approximately logarithmically 
along the real line. This means that, strictly speaking, all sets representable by 
the computer are closed, compact and totally disconnected. Continuity cannot be 
sensibly defined for functions on a totally disconnected domain, and many functions, 
such as y — |x, are, surprisingly, not one-to-one. In this example, two numbers 
adjacent in the representation, when multiplied by |, may both produce the same 
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Figure 2. Schematic depiction of a geodesic crossing from the chart {Ua,4'a) to the chart {Ug,4>0)- 

floating-point result, as the difference between the true results is less than the 

precision of the discrete representation. 

Thus, although a computer cannot directly represent a manifold, atlas or chart 
as defined above, it can produce similar objects, which we distinguish from the 
abstract entities by use of a Sans Serif typeface. 

We define an Atlas, the numerical representation of an atlas {(Z^aj'^a)}) as a 
collection of Charts, the numerical representations of charts {Ua,(l>a)- We define 
the Chart representing {Ua,(l^a) as: 

1. The numerical representation of the set (j)aipia)- 

2. The set of numerical representations of functions mapping from the Chart to 
other Charts: {(f>i:i o (f>-^ : (f>a{Uar]Up) —> (f)fj{Ua nUp)}. 

We cannot define a Manifold as an Atlas combined with the set of equivalence classes 
of points X e (j}a{^a) Under the mappings (^/j o as these mappings are not, in 
general, one-to-one, due to the discrete representation employed by the computer. 
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It is possible that 0q, o (p'j^'^ {(j) o 4'a'^{x)) ^ x, though the difference should be 
"small", that is, on the order of the local resolution e of the discrete representation. 
In general, there appears to be no sensible definition for a Manifold, and we do not 
adopt one. 

Although we cannot construct an authoritative set A4 using the above proce- 
dure, we can still produce a useful definition of a Point as the numerical represen- 
tation of p e A4. First define a Coordinate^ by: 

1. The numerical representation (Chart) of a chart a (abbreviated notation for 

2. The numerical representation of a point x £ (j)a{l^a)- 

Note that the numerical representation of x is simply an n-tuple of real numbers, 
so it is necessary to additionally specify the chart in order to give those numbers 
the context of a mapping. We can identify a Coordinate {a,x) with a point of the 
manifold, 4>~^(x) £ Ua C Ai. We do not, however, have a numerical representation 
of the function (f>~^, so to define a Point representing p £ Ua use a set of 
Coordinates: 

1. The Coordinate {a,(t>a{p) — x). 

2. The Coordinates {(/3, 0^ o 0^i((/)„(p))) : p^a}. 

Note that these are not the Coordinates (/3, (f)i3{p)), though the difference should be 
"small". 

Thus, although we can define a Point, its definition is tied to its coordinates 
under a particular chart. If we were to produce Coordinates (a, (paip)) and (/3, 4'fj{p)) 
for the same point p £ Ua C)Ui3 C Ai, the Points manufactured from them would 
not, in general, consist of precisely the same Coordinates (condition 2 above), and 
thus would not be equal so far as the computer is concerned. We cannot produce a 
chart-independent representation of a point; we can only produce a chart-dependent 
Point that is "approximately" chart-independent, in the sense that the difference 
between the representations arising from different charts is "small" . 

There is an obvious similarity between the Atlas and Point; both are collections 
of different numerical representations of the same abstract object. This abstract 
layer identifying different numerical representations imparts (approximate) chart- 
independence to numerical operations performed within the structure. Essentially, 
in the midst of performing a computation the computer can change charts as re- 
quired by mapping the pertinent parameters to another Chart. 

An example we will follow throughout this paper is the computation of a discrete 
approximation to a geodesic. A geodesic 7 : R — > may, for some interval around 
r e R, lie in the domain Ua of a chart a. It can be represented numerically 
by the series of Coordinates (and thus by the series of Points generated from the 
Coordinates) (a, (j)a{^{T))), {a, (j)ai'yiT + S))), (a, 4>ai7iT + 26))), . . . for some step- 
size 6 > 0, but it is possible that for some k £ N, ^{t + (fc -|- 1)6) ^ Ua- In 
this case, for 5 chosen to be sufficiently small, there exists some Chart (3 such that 

"^Internally named a Node by GRworkbench for historical reasons. 
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7(t + (fc + 1)5) e Uf3, and j{t + kS) G Ua r]Up. Without the abstract identification 
of the Coordinate (a, (j)a{j{T + k6))) with the Coordinate (/3, 0/3o0~^((^Q,(7(r + fc(5)))) 
via the Point representing the point 7(t + fc(5), there would be no way in which the 
next Point, representing the point 7(t + {k + 1)6), could be computed. With the 
identification, however, we can continue to compute the numerical representation 
from algorithms operating on Chart f3. This procedure is depicted schematically in 



The interest of users in entities, such as geodesies, that are structured sets of 
Points, motivates us to define the Object. An Object is a set of Points. The Points 
constituting an Object will have Coordinates. We define a Segment for a certain 
Object and Chart as the image of the set represented by the Object on that Chart. 
As such, a Segment is a set of Coordinates. 

To use the differential geometric framework to represent a space-time, we must 
provide one additional component. We extend the definition of a Chart to encompass 
the provision of a metric g — a (0, 2) symmetric tensor field on M. We require that 
each Chart provide a function returning the tensor components gij\x with respect 
to the chart's coordinate basis {d/dx"'}, for any coordinate x G ^q(Z^q). 

3 Implementation 

The numerical differential geometric system defined above may be naturally ex- 
pressed as a collection of classes in the C-|--fQ programming language. A class is 
itself a collection of data and related functions. The notation A::B indicates that B 
is a member of class A. 

The Atlas class stores chart-independent data and a list of the Charts {{Ua, (f>a)} 
that comprise it. In the case of the Atlas representing the Schwarzschild space- 
time, the Atlas stores the Schwarzschild mass M that defines certain properties of 
the Charts, such as mandating that the radius xi > 2M for exterior charts. 

The Chart representing {lAa,4>a) must provide 

1. The set (j)a{Ua) C R". 

2. The functions {(jyp o : (jy^ipla CiUp) (pfsiUa CiUfs)}. 

3. The metric tensor gij\x, for x £ (pa^a)- 

For maximum fiexibility, we require the user to define the set (j)a (Ma ) in terms of a 
Boolean function on M", that is, a function Chart::lnterior : R" {true, false}. 



This formalism allows the user to plug in any algorithm to define the domain of the 
Chart. 

For an exterior spherical polar chart of the Schwarzschild space-time, the 
Chart::lnterior function is defined as 



Figure 



I- 




(1) 



Chart::lnterior(a;) = 



true if xi > 2M and < a;2 < tt and < X3 < 27r 
false otherwise 



(2) 
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Maps to other charts, 0/3 o are implemented as M" M" functions. An 
obvious disadvantage of the need to define maps (j)0o^~'^ rather than simply maps 
(j)a is that, in general, an Atlas comprising m Charts requires up to m,{m — 1) 
maps rather than just the m A4 ^ IR." maps. This number can be 
reduced by permitting the implicit definition of maps, by application of several 
maps between other charts, such as the map 0/3 o 0"-"^ = (0^ o 0~"^) o (0^ o : 

Two spherical polar charts are required to cover the entire exterior [xi > 2M) 
region of the Schwarzschild space-time. Choosing the second spherical polar chart 



scott-grwb: submitted to World Scientific on February 7, 2008 6 



so that its polar axis corresponds to the coordinates^ X2 = 7r/2 and X3 = ■n/2,'i'K/2 
allows us to define the map between the charts as: 



x'l — Xi 

x'2 — arg(sinx2 sinxa + \/ sin^ X2 sin^ X3 
X3 = arg(— sina;2 C0SX3 + ^/—l cosa;2) 



(3) 



This map is its own inverse. 

The numerical representation of a differential geometric structure cannot know, 
in advance, if that representation is complete. That is, the implementation cannot 
determine if the Atlas provided covers the entire manifold A4, as the implementation 
has no a priori knowledge of the manifold M itself. Neither can the implementation 
determine if the E" R" maps supplied by the Charts are comprehensive, as the 
relationships between Charts are defined solely in terms of these maps. This allows 
the user freedom to implement systems for which the Atlas provided does not cover 
the manifold, and for which the M" M" maps supplied by the Charts are not 
comprehensive. 

This is not cause for concern. When the implementation cannot successfully 
continue correct computation due to the lack of a Chart or map, the algorithm halts. 
If the user wishes to perform the computation, they must add a new component to 
the incomplete system. In many circumstances, though, the user may be concerned 
only with a portion of the system. For example, to study orbits in a Schwarzschild 
space-time, it is sufficient to define only exterior charts. To require the user to define 
portions of the space-time which they have no intention of utilizing — in this case, 
the interior Schwarzschild space-time — would be to unnecessarily inconvenience the 
user. 

To complete the definition of the Chart, the user must supply a function return- 
ing the metric tensor components gij\xi for x € (j)a{Ua). For any exterior spherical 
polar chart of the Schwarzschild space-time, the components are defined as: 

500U = -1 + ^ .9oiU = 502U = ffo3U = 

gioU = o giiU = (i-^)~' 5i2U = o 5i3U = o 

g2o|x = g2i\x = 522U = 523 U = 

ffsoU^O g3i|x = 53211 = 533U = (a;i)^sin^X2 

In the current implementation. Atlases and Charts are defined by the user by 
modifying trivial C++ code in supplied "template" file^. Only the most basic 
programming skills are required of the user. Compared to an obvious alternative, 
namely, using a scripting language, this method has the advantage of producing 
efficient, pre-compiled code, but the disadvantage of relinking the executable when 
space-time definitions are modified. Such modifications are, however, comparatively 
rare for common usage. 

The remaining classes. Objects, Segments, Points and Coordinates are imple- 
mented in a straightforward way. Coordinates contain an n-tuple of floating point 
values and an indirect reference (via their containing Segment) to their Charts. Any 



^ We use the C-l — h convention that indices begin from 0. 
'^Not to be confused with the CH — h template keyword. 
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Point owns a list of its Coordinate representations on all applicable Charts. A Seg- 
ment similarly maintains a list of the Coordinates of the Points of a particular Object 
on a particular Chart. An Object maintains a list of its Points, and Segments on 
various Charts, and a Chart maintains a list of Segments on it from various Ob- 
jects. Segments and Coordinates thus belong to two lists, which are represented 
perpendicularly to one another in Figure ^. 

The above components are genuinely sufficient to define a space-time. Where 
GRworkbench requires certain information to perform a task, such as the geodesic 
equation to compute a geodesic, it extracts that information numerically using the 
relevant definition. For example, the geodesic equation is given by 

d^x-^ _ dxUx- 
dr2 - ^ dr dr ' 

where the Christoffel symbol F is given by 

We may compute as the inverse matrix of gab, and compute the partial deriva- 
tives -^gab by numerical differentiation of gab- Thus, GRworkbench can compute 
the geodesic equation directly from the user-supplied metric for any given space- 
time. 



4 Utility 

Using the definition for a Schwarzschild Atlas given in Section 3, consisting of two 
exterior spherical polar Charts, we can begin to trace geodesies, such as those of 
Figure 0, on the space-time without deriving the geodesic equations. 

A well-known fact about the Schwarzschild space-time is that there exists a 
circular null orbit at xi = 3M. Although this fact is not immediately apparent 
analytically, it can be reproduced using GRworkbench in an experimental mode. 

We release geodesies from xq — 0,X2 — 7r/2,X3 = tt/2 for varying values of 
xi > 2M, and mandate that ^ > 0, ^ = ^ = 0, ^ > and the geodesic 
be null. GRworkbench then uses the metric to produce an initial null tangent 
vector satisfying these condition^. For any xi, the geodesic will either escape to 
infinity, or fall into the event horizon. Using these conditions, we can iterate down 
on the value of xi that divides these two types of behaviour. Figure ^ shows the 
geodesies traced in this process. As expected, the value below which geodesies fall 
into the event horizon, and above which geodesies escape to infinity, is given by 
xi = 1.5 = 3 X 0.5 = 3M to high accuracy. The computations were performed 
interactively in a few minutes on an inexpensive notebook computer. 

While this is a trivial example, recall that this functionality is available purely 
from a minimal definition of the space-time — and as such, is available for any space- 
time that can be so defined. Additionally, there are many algorithms implemented, 

''GRworkbench does this by breaking the given vector into a purely time-like and purely space-like 
component, and re-scaling these parts so that they sum to a null vector. 
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Figure 4. Null geodesies in the Schwcirzschild spcLce-time, iterating in on the r = 1.5 = 3M circular 
null orbit. 



beyond the featured example of geodesic tracing, that operate; witli similar flexibil- 
ity; from the production of null cones to the examination of causal connections, to 
the computation of proper distances between points on a space-like hypersurface. 
GRworkbench's firm basis in differential geometry makes it a truly general tool, and 
encourages a new experimental approach to problems in General Relativity — one 
of trial and observation. 



5 Conclusion 

GRworkbench successfully implements a numerical analogue of a standard differ- 
ential geometric system, mirroring a manifold and atlas of charts with the CH — h 
classes Atlas and Chart. Provision of a metric on a Chart completes the definition 
of a space-time in General Relativity. 

From this minimal definition, complex algorithms, such as geodesic tracing, 
may be "bootstrapped" through layers of numerical operations, to produce useful 
results, in novel ways, with minimal effort. As such, GRworkbench permits, and 
encourages, an "experimental" approach to problem-solving in General Relativity. 
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